#**********************************************************************************
# ccamSens.r
# This file contains the code necessary to plot sensitivities overlaid on top
# of one another by sensitivity group which is set in the file 'SensitivityGroup'
# in each Scenario folder.  Scenarios with the same sensitivity group will be plotted
# together.  The 'SensitivityGroup' file should have a single number in it, nothing more.
# The Base model will be SensitivityGroup=0 and will be plotted against all other
# Sensitivity groups.
#
# Assumes the opList has been built and contains the SensitivityGroup element:
# opList[[scenario]][[4]]$SensitivityGroup
#
# Author            : Chris Grandin
# Development Date  : December 2011 - September 2012
#
#**********************************************************************************

fig.base.vs.sens <- function(sensitivityGroup=1,whichPlot="biomass",ylimit=6,useMaxYlim=T,lty=2,lwd=2,pch=20,offset=0.3,opacity="20"){
  # plots Spawning stock biomiass, depletion, and recruitment for a given sensitivity group
  #  - plot the MCMC posterior data for each, with confidence limits
  #  - offset is the number of years to offset each vertical bar from each other.
  #  - opacity is a two digit string from 00-99
  # whichPlot can be:
  # 1. "biomass"
  # 2. "depletion"
  # 3. "recruits"

  try(dev.off(), silent=T)
  op	<- par(no.readonly=T)
  base <- 0
  color <- 1
  colors <- color

  for(scenario in 1:length(opList)){
    if(opList[[scenario]][[4]]$SensitivityGroup == 0){
      base <- scenario
    }
  }
  plotR <- T
  baseRep <- opList[[base]][[4]]
  runNames <- strsplit(opList[[base]][[1]],"/")[[1]][3] # Gets the run's folder name out of the reletive path
  if(whichPlot == "recruits"){
    if(opList[[base]][[7]]){ # if MCMC results are loaded for base
      mc <- baseRep$mc.rt
      mc.rt <- as.data.frame(window(mcmc(mc),start=Burn,thin=Thin))
      rt <- apply(mc.rt,2,quantile,probs=c(0.025,0.5,0.975))
      if(useMaxYlim){
        yUpperLimit <- max(rt)
      }else{
        yUpperLimit <- ylimit
      }
      xp <- plot(ryr,
                 rt[2,],
                 type="p",
                 pch=20,
                 col=color,
                 xlim=c(min(ryr),max(ryr)),
                 ylim=c(0,yUpperLimit),
                 xlab="Year",
                 ylab="Age-1 recruits (billions)",
                 main="Age-1 recruits",
                 las=1)
      arrows(ryr, rt[1, ],ryr,rt[3,],code=3,angle=90,length=0.01,col=color)
      abline(h=median(as.matrix(mc.rt)),col=2,lty=2)
      abline(h=mean(as.matrix(mc.rt)),col=3,lty=2)
      currOffset <- offset
      for(scenario in 1:length(opList)){
        if(scenario != base){
          if(opList[[scenario]][[4]]$SensitivityGroup == sensitivityGroup){
            if(opList[[scenario]][[7]]){ # if MCMC results are loaded
              mc <- opList[[scenario]][[4]]$mc.rt
              mc.rt <- as.data.frame(window(mcmc(mc),start=Burn,thin=Thin))
              rt <- apply(mc.rt,2,quantile,probs=c(0.025,0.5,0.975))
              color <- color + 1
              colors <- c(colors,color)
              par(new=T)
              plot(ryr+currOffset,
                   rt[2,],
                   pch=20,
                   xlim=c(min(ryr),max(ryr)),
                   ylim=c(0,yUpperLimit),
                   col=color,
                   xlab="",
                   ylab="",
                   las=1,
                   axes=F)
              arrows(ryr+currOffset, rt[1, ],ryr+currOffset,rt[3,],code=3,angle=90,length=0.01,col=color)
              par(new=F)
              currOffset <- currOffset + offset
              runNames <- c(runNames,strsplit(opList[[scenario]][[1]],"/")[[1]][3])
            }else{
              cat(paste("Error plotting 'Recruitment Sensitivity' figure.  There are no MCMC outputs loaded for scenario",scenario,", so it is not shown on the plot.\n\n"))
            }
          }
        }
      }
      lty <- c(rep(1,(length(runNames))),2,2)
      runNames <- c(runNames,"long-term median","long-term mean")
      legend("topright",
             runNames,
             lty=lty,
             col=c(colors,2,3),
             bty="n",
             lwd=2)
      filename <- paste("SensitivityGroup_",sensitivityGroup,"_Recruits",sep="")
      saveFig(filename)
    }else{
      cat(paste("Error plotting 'Recruitment Sensitivity' figure.  There are no MCMC outputs loaded for the base scenario (",base,")\n\n"))
    }
  }else if(whichPlot == "biomass"){
    if(opList[[base]][[7]]){ # if MCMC results are loaded for base
      mcbo <- baseRep$mc$bo
      post.bo <- as.data.frame(window(mcmc(mcbo),start=Burn,thin=Thin))
      boci <- apply(post.bo,2,quantile,probs=c(0.025,0.5,0.975))

      mcbt <- baseRep$mc.sbt
      post.bt <- as.data.frame(window(mcmc(mcbt),start=Burn,thin=Thin))
      btci <- apply(post.bt,2,quantile,probs=c(0.025,0.5,0.975))

      if(useMaxYlim){
        yUpperLimit <- max(btci)
      }else{
        yUpperLimit <- ylimit
      }
      matplot(baseRep$yrs,
              t(btci),
              type="l",
              col=color,
              lty=c(2,1,2),
              lwd=2,
              ylim=c(0,yUpperLimit),
              xlab="Year",
              ylab="Spawning biomass",
              main="Spawning biomass",
              las=1)

      # Shade the confidence interval
      xx <- c(baseRep$yrs,rev(baseRep$yrs))
      yy <- c(btci[1,],rev(btci[3,]))
      shade <- getShade(color,opacity)
      polygon(xx,yy,density=NA,col=shade)
      # End shade the confidence interval

      points(baseRep$yrs[1]-0.8,boci[2],col=color,pch=1)
      arrows(baseRep$yrs[1]-0.8,boci[1],baseRep$yrs[1]-0.8,boci[3],col=color, code=0, lwd=1.5)
      currOffset <- offset
      for(scenario in 1:length(opList)){
        if(scenario != base){
          if(opList[[scenario]][[4]]$SensitivityGroup == sensitivityGroup){
            if(opList[[scenario]][[7]]){ # if MCMC results are loaded
              mcbo <- opList[[scenario]][[4]]$mc$bo
              post.bo <- as.data.frame(window(mcmc(mcbo),start=Burn,thin=Thin))
              boci <- apply(post.bo,2,quantile,probs=c(0.025,0.5,0.975))

              mcbt <- opList[[scenario]][[4]]$mc.sbt
              post.bt <- as.data.frame(window(mcmc(mcbt),start=Burn,thin=Thin))
              btci <- apply(post.bt,2,quantile,probs=c(0.025,0.5,0.975))

              color <- color + 1
              colors <- c(colors,color)
              par(new=T)
              matplot(opList[[scenario]][[4]]$yrs,
                      t(btci),
                      type="l",
                      col=color,
                      lty=c(2,1,2),
                      lwd=2,
                      ylim=c(0,yUpperLimit),
                      xlab="",
                      ylab="",
                      las=1)

              # Shade the confidence interval
              xx <- c(opList[[scenario]][[4]]$yrs,rev(opList[[scenario]][[4]]$yrs))
              yy <- c(btci[1,],rev(btci[3,]))
              shade <- getShade(color,opacity)
              polygon(xx,yy,density=NA,col=shade)
              # End shade the confidence interval

              par(new=F)
              points(opList[[scenario]][[4]]$yrs[1]-0.8+currOffset,boci[2],col=color,pch=1)
              arrows(opList[[scenario]][[4]]$yrs[1]-0.8+currOffset,boci[1],A$yrs[1]-0.8+currOffset,boci[3],col=color, code=0, lwd=1.5)
              currOffset <- currOffset + offset
              runNames <- c(runNames,strsplit(opList[[scenario]][[1]],"/")[[1]][3])
             }else{
              cat(paste("Error plotting 'Spawning Biomass Sensitivity' figure.  There are no MCMC outputs loaded for scenario",scenario,", so it is not shown on the plot.\n\n"))
            }
          }
        }
      }
      lty <- c(rep(1,(length(runNames))))
      legend("topright",runNames,lty=lty,col=colors,bty="n", lwd=2)

      filename <- paste("SensitivityGroup_",sensitivityGroup,"_Biomass",sep="")
      saveFig(filename)
    }else{
      cat(paste("Error plotting 'Spawning Biomass Sensitivity' figure.  There are no MCMC outputs loaded for the base scenario (",base,")\n\n"))
    }
  }else if(whichPlot == "depletion"){
    if(opList[[base]][[7]]){ # if MCMC results are loaded for base
      mcdt <- baseRep$mc.sbdepletion
      post.dt  <- as.data.frame(window(mcmc(mcdt),start=Burn,thin=Thin))
      dtci <- apply(post.dt,2,quantile,probs=c(0.025,0.5,0.975))
      if(useMaxYlim){
        yUpperLimit <- max(dtci)
      }else{
        yUpperLimit <- ylimit
      }
      matplot(baseRep$yrs,
              t(dtci),
              type="l",
              col=color,
              lty=c(2,1,2),
              lwd=2,
              ylim=c(0,yUpperLimit),
              xlab="Year",
              ylab="Spawning depletion",
              main="Depletion")

      # Shade the confidence interval
      xx <- c(baseRep$yrs,rev(baseRep$yrs))
      yy <- c(dtci[1,],rev(dtci[3,]))
      shade <- getShade(color,opacity)
      polygon(xx,yy,density=NA,col=shade)
      # End shade the confidence interval

      abline(h=0.40, lwd=mtLineWidth, col=mtLineColor, lty=mtLineType)
      for(scenario in 1:length(opList)){
        if(scenario != base){
          if(opList[[scenario]][[4]]$SensitivityGroup == sensitivityGroup){
            if(opList[[scenario]][[7]]){ # if MCMC results are loaded
              mcdt <- opList[[scenario]][[4]]$mc.sbdepletion
              post.dt  <- as.data.frame(window(mcmc(mcdt),start=Burn,thin=Thin))
              dtci <- apply(post.dt,2,quantile,probs=c(0.025,0.5,0.975))
              color <- color + 1
              colors <- c(colors,color)
              par(new=T)
              matplot(opList[[scenario]][[4]]$yrs,
                      t(dtci),
                      type="l",
                      col=color,
                      lty=c(2,1,2),
                      lwd=2,
                      ylim=c(0,yUpperLimit),
                      xlab="",
                      ylab="")

              # Shade the confidence interval
              xx <- c(opList[[scenario]][[4]]$yrs,rev(opList[[scenario]][[4]]$yrs))
              yy <- c(dtci[1,],rev(dtci[3,]))
              shade <- getShade(color,opacity)
              polygon(xx,yy,density=NA,col=shade)
              # End shade the confidence interval

              par(new=F)
              runNames <- c(runNames,strsplit(opList[[scenario]][[1]],"/")[[1]][3])
            }else{
              cat(paste("Error plotting 'Depletion Sensitivity' figure.  There are no MCMC outputs loaded for scenario",scenario,", so it is not shown on the plot.\n\n"))
            }
          }
        }
      }
      lty <- c(rep(1,(length(runNames))),mtLineType)
      lwd <- c(rep(2,(length(runNames))),mtLineWidth)
      colors <- c(colors,mtLineColor)
      runNames <- c(runNames,"Management Target")
      legend("topright",runNames,lty=lty,col=colors,bty="n", lwd=lwd)
      filename <- paste("SensitivityGroup_",sensitivityGroup,"_Depletion",sep="")
      saveFig(filename)
    }else{
      cat(paste("Error plotting 'Depletion Sensitivity' figure.  There are no MCMC outputs loaded for the base scenario (",base,")\n\n"))
    }
  }
  par(op)
}
